/*
 * Simulation.cpp
 *
 *  Created on: Jan 30, 2012
 *      Author: falko
 */


#include "Simulation.h"

#include <algorithm>
#include <ctype.h>
#include <fstream>

namespace falko
{


void ParticleFactory::Generate(particle * p, particle_type type)
{
	switch (type)
	{
	case ELECTRON: electron(p); break;
	case PROTON:   proton(p);   break;
	case H_PLUS:   h_plus(p);   break;
	default:       clear(p);    break;
	}
}

void ParticleFactory::clear(particle * p)
{
	p->type_ = DEFAULT;
	p->id_ = -1;
	p->mass_ = 0.0;
	p->abs_charge_ = 0.0;
	p->charge_ = 0;
	p->position_ = 0.0;
	p->velocity_ = 0.0;
	p->acceleration_ = 0.0;
}

void ParticleFactory::electron(particle * p)
{
	p->type_ = ELECTRON;
	p->id_ = nextID();
	p->mass_ = "9.109e-31";
	p->abs_charge_ = "1.60217646e-19";
	p->charge_ = -1;
	p->position_ = 0.0;
	p->velocity_ = 0.0;
	p->acceleration_ = 0.0;
}

void ParticleFactory::proton(particle * p)
{
	p->type_ = PROTON;
	p->id_ = nextID();
	p->mass_ = "1.67262158e-27";
	p->abs_charge_ = "1.60217646e-19";
	p->charge_ = 1;
	p->position_ = 0.0;
	p->velocity_ = 0.0;
	p->acceleration_ = 0.0;
}

void ParticleFactory::h_plus(particle * p)
{
	p->type_ = H_PLUS;
	p->id_ = nextID();
	p->mass_ = "3.34524316e-27";
	p->abs_charge_ = "3.2043529200000001e-19";
	p->charge_ = 1;
	p->position_ = 0.0;
	p->velocity_ = 0.0;
	p->acceleration_ = 0.0;
}






SimulationFileWriter::~SimulationFileWriter()
{
	for (WriterMap::const_iterator it = writers_map_.begin(); it != writers_map_.end(); ++it)
	{
		delete it->second;
	}
}

void SimulationFileWriter::write (particle * p, hpf_t const &time)
{
	if (!writers_map_.count(p))
	{
		writeNewParticle(p);
	}
	writers_map_[p]->vel_file_ << time << "\t" << p->velocity_ << "\n";
	writers_map_[p]->pos_file_ << time << "\t" << p->position_ << "\n";
}

void SimulationFileWriter::writeNewParticle(particle * p)
{
	std::stringstream pos;
	pos << simulation_name_ << "_" << p->id_ << "_position"<< ".dat";
	std::stringstream vel;
	vel << simulation_name_ << "_" << p->id_ << "_velocity"<< ".dat";

	particle_writers * writer = writers_map_[p] =
	            new particle_writers(pos.str().c_str(), vel.str().c_str());

	if (!(writer->pos_file_ && writer->vel_file_))
	{
		std::cerr << "bad file...\n";
		exit (1);
	}

	writer->pos_file_ << "#\t" << pos.str().c_str() << "\n";
	writer->vel_file_ << "#\t" << vel.str().c_str() << "\n";
	writeParticleData(p, writer->pos_file_);
	writeParticleData(p, writer->vel_file_);
	writer->pos_file_ << "time\tx-axis\ty-axis\tz-axis\n";
	writer->vel_file_ << "time\tx-axis\ty-axis\tz-axis\n";
}

void SimulationFileWriter::writeParticleData(particle * p, HpFileWriter &w)
{
	w << "#\t" << p->mass_ << "\n";
	w << "#\t" << p->abs_charge_ * p->charge_ << "\n";
}



ParticleSource::ParticleSource(particle_type t, hpf_t v,
								uint max, hpv3d_t c, hpv3d_t dir)
: type_(t), particles_(max), voltage_(v),
  max_emissions_(max), center_(c), emision_dir_(dir)
{
	emision_dir_.normalize();
}

ParticleSource::~ParticleSource()
{
	for (int i = particles_.size() - 1; i >= 0; --i)
	{
		delete particles_[i];
	}
}

particle * ParticleSource::getParticle()
{
	particle * temp = new particle(type_);
	hpf_t velocity_mag = (( temp->abs_charge_ * 2.0 * voltage_) / temp->mass_).sqrt();
	temp->velocity_ = emision_dir_ * velocity_mag ;
	temp->position_ = center_;
	particles_.push_back(temp);
	max_emissions_--;
	return temp;
}







Simulation::~Simulation()
{
	if (ion_source_ != NULL)
		delete ion_source_;
}

void Simulation::run(const char * file)
{
	parseFile(file);

	switch (type_)
	{
	case GENERATE_VECTOR_FIELD:
		generateVectorField();
		break;
	case PARTICLE_TRACE:
		particleTrace();
		break;
	default:
		break;
	}

}

void Simulation::particleTrace()
{
	if (ion_source_ == NULL)
		return;

	HpFileWriter sout;
	std::cout << "running:     " << simulation_name_ << std::endl;
	std::cout << "Iterations:  " << iterations_ << std::endl;
	sout      << "time step:   " << delta_time_ << "\n";

	for (int j = iterations_; j > 0; --j)
	{
		step();
	}
}

void Simulation::generateVectorField()
{
	std::stringstream ss;
	ss << simulation_name_ << "_mag_field.dat";
	{
		HpFileWriter magFile(ss.str().c_str());
		if(!(magFile))
		{
			std::cerr << "could not open file" << std::endl;
			return;
		}

	//	TODO: center around all coils ... remove magic numbers
		hpv3d_t point(0.0, 0.0, 0.0);
		hpv3d_t fieldLine(0.0, 0.0, 0.0);
		hpf_t step(0.0);
		const double initCoord = -.5;
		const double pointCount = 10;

		point = initCoord;
		step = (fabs(initCoord)*2) / pointCount;

		magFile << "#Px\tPy\tPz\tBx\tBy\tBz" << "\n";
		for (uint z = 0; z < pointCount; ++z)
		{
			for (uint y = 0; y < pointCount; ++y)
			{
				for (uint x = 0; x < pointCount; ++x)
				{
					mag_field_.getField(point, fieldLine);
					magFile << point << "\t" << fieldLine << "\n";
					point.x_ += step;
				}
				point.x_ = initCoord;
				point.y_ += step;
			}
			point.y_ = initCoord;
			point.z_ += step;
		}
	}
	plotVectorField(ss.str().c_str());
}

void Simulation::step()
{
	//inject particles
	if ((*ion_source_))
		particles_.push_back(ion_source_->getParticle());

	//update particles
	hpv3d_t b_field(0.0,0.0,0.0);
	hpv3d_t k1;
	hpv3d_t k2;
	hpv3d_t k3;
	hpv3d_t k4;
	hpv3d_t temp_velocity;

	for (int i = particles_.size() - 1; i >= 0; --i)
	{
		particle* pa = particles_[i];
		pa->position_ += pa->velocity_ * delta_time_;
		temp_velocity = pa->velocity_;
		mag_field_.getField(pa->position_, b_field);

		calculateAcceleration(pa, b_field, k1);
		k1 *= delta_time_;

		pa->velocity_ = temp_velocity + (k1 / 2.0);
		calculateAcceleration(pa, b_field, k2);
		k2 *= delta_time_;

		pa->velocity_ = temp_velocity + (k2 / 2.0);
		calculateAcceleration(pa, b_field, k3);
		k3 *= delta_time_;

		pa->velocity_ = temp_velocity + k3;
		calculateAcceleration(pa, b_field, k4);
		k4 *= delta_time_;

		pa->velocity_ = temp_velocity + ((k1 + k2 + k3 + k4) / 6.0);

//		calculateAcceleration(pa, b_field, k1);
//		pa->velocity_ += k1 * delta_time_;

		pa->velocity_.normalize();
		pa->velocity_ *= temp_velocity.magnitude();

		file_writter_.write(pa, time_);
	}
	time_ += delta_time_;
}

void Simulation::calculateAcceleration(particle const * pa,
						   hpv3d_t const &b_field,
						   hpv3d_t &acceleration)
{
	hpv3d_t b_field_unit;
	b_field_unit = b_field;
	b_field_unit.normalize();
	acceleration= (
			(pa->velocity_ - b_field_unit *
			(pa->velocity_ * b_field_unit)) % b_field)
			* ((pa->abs_charge_  / pa->mass_) * pa->charge_);
}


#define IS_ALPHA(c)    (((c >= 65) && (c <= 172)) && ((c < 91 || c > 96)))
#define IS_NUMERIC(c)  (((c >= 48) && (c <= 57)) || (c == '.') || (c == '-') || (c == 'E'))
//int not_alpha_numeric(char c) { return !(IS_ALPHA(c) || IS_NUMERIC(c)); }
int not_alpha(char c) { return !IS_ALPHA(c); }
std::string extract(std::stringstream &ss, int delim)
{
	std::stringstream iss;
	while (ss.peek() != delim && !ss.fail() )
	{
		iss << static_cast<char>(ss.get());
	}
	return iss.str();
}

hpf_t get_hpf(std::stringstream &ss)
{
	for (;!ss.fail() && !IS_NUMERIC(ss.peek()); ss.get());
	std::stringstream num;
	while (!ss.fail() && IS_NUMERIC(ss.peek()))
	{
		num << static_cast<char>(ss.get());
	}
	return hpf_t(num.str().c_str());
}

hpv3d_t get_hpv3d(std::stringstream &ss)
{
	hpv3d_t num;
	num.x_ = get_hpf(ss);
	num.y_ = get_hpf(ss);
	num.z_ = get_hpf(ss);
	return hpv3d_t(num);
}

void* get_ini_line(std::ifstream &ifs, std::string &line)
{
	void* ret = 0;
	if ((ret = getline(ifs, line)))
	{
		while ((line.empty() || line.at(0) == '#') && (ret = getline(ifs, line)))
			;
		if (ret)
		{
			line.erase(remove_if(line.begin(), line.end(), ::isspace), line.end());
			std::transform(line.begin(), line.end(),line.begin(), ::toupper);
		}
	}
	return ret;
}

void Simulation::parseFile(const char * file)
{
	std::ifstream ifs(file);
	if (!ifs.is_open())
	{
		std::cerr << "could not open file" << std::endl;
		return;
	}

	std::string line;
	std::vector<std::vector<std::string>> component_list;
	while (get_ini_line(ifs, line))
	{
		if (!IS_ALPHA(line.at(0)))
		{
			std::vector<std::string> component;
			line.erase(remove_if(line.begin(), line.end(), not_alpha), line.end());
			std::string tag = line;
			component.push_back(tag);
			while (get_ini_line(ifs, line) && IS_ALPHA(line.at(0)))
			{
				component.push_back(line);
			}
			line.erase(remove_if(line.begin(), line.end(), not_alpha), line.end());
			if (tag != line)
			{
				std::cout << tag << std::endl;
				std::cout << line << std::endl;
				std::cerr << "File Error: Tag Match" << std::endl;
				ifs.close();
				return;
			}
			component_list.push_back(component);
		}
		else
		{
			std::cerr << "File Error: Missing Tag" << std::endl;
			ifs.close();
			return;
		}
	}
	ifs.close();

	for (uint i = 0; i < component_list.size(); ++i)
	{
		if(component_list[i].at(0) == "SIMULATION")
			initSimulation(component_list[i]);
		else
			if (component_list[i].at(0) == "COIL")
				addCoil(component_list[i]);
			else
				if (component_list[i].at(0) == "PARTICLESOURCE")
					initIonSource(component_list[i]);

	}
}

void Simulation::initSimulation(std::vector<std::string> &init_data)
{
	time_ = 0.0;
	for (uint i = 1; i < init_data.size(); ++i)
	{
		std::stringstream ss(init_data[i]);
		std::string var = extract(ss, '=');
		ss.get();

		if (var == "NAME")
		{
			std::string s;
			ss >> s;
			simulation_name_ = s;
			file_writter_.setName(s);
		}
		else if (var == "TYPE")
		{
			std::string var;
			ss >> var;
			if (var == "GENERATEVECTORFIELD")
			{
				type_ = GENERATE_VECTOR_FIELD;
			}
			else if (var == "PARTICLETRACE")
			{
				type_ = PARTICLE_TRACE;
			}
		}
		else if (var == "ITERATIONS")
		{
			iterations_ = static_cast<uint>(get_hpf(ss).toDouble());
		}
		else if (var == "DELTATIME")
		{
			delta_time_ = get_hpf(ss);
		}
	}
}

void Simulation::addCoil(std::vector<std::string> &init_data)
{
	uint    loops;
	hpf_t   wire_radius;
	hpf_t   length;
	hpf_t   inner_radius;
	hpf_t   current;
	hpv3d_t center;

	HpFileWriter sout;

	for (uint i = 1; i < init_data.size(); ++i)
	{
		std::stringstream ss(init_data[i]);
		std::string var = extract(ss, '=');
		ss.get();

		if (var == "LOOPS")
		{
			loops = static_cast<uint>(get_hpf(ss).toDouble());
		}
		else if (var == "WIRERADIUS")
		{
			wire_radius = get_hpf(ss);
		}
		else if (var == "LENGTH")
		{
			length = get_hpf(ss);
		}
		else if (var == "INNERRADIUS")
		{
			inner_radius = get_hpf(ss);
		}
		else if (var == "CURRENT")
		{
			current = get_hpf(ss);
		}
		else if (var == "CENTER")
		{
			center = get_hpv3d(ss);
		}
	}

	PathFactory::MakeCoil(&mag_field_,loops, wire_radius,
			length, inner_radius, current, center);
}

void Simulation::initIonSource(std::vector<std::string> &init_data)
{
	if (ion_source_ != NULL)
		return;

	particle_type type;
	hpf_t voltage;
	uint max_emissions = 0;
	hpv3d_t center;
	hpv3d_t dir;

	for (uint i = 1; i < init_data.size(); ++i)
	{
		std::stringstream ss(init_data[i]);
		std::string var = extract(ss, '=');
		ss.get();

		if (var == "TYPE")
		{
			std::string var;
			ss >> var;
			if (var == "ELECTRON")
			{
				type = ELECTRON;
			}
			else if (var == "PROTON")
			{
				type = PROTON;
			}
			else if (var == "H_PLUS")
			{
				type = H_PLUS;
			}
		}
		else if (var == "VOLTAGE")
		{
			voltage = get_hpf(ss);
		}
		else if (var == "MAXPARTICLES")
		{
			max_emissions = static_cast<uint>(get_hpf(ss).toDouble());
		}
		else if (var == "CENTER")
		{
			center = get_hpv3d(ss);
		}
		else if (var == "EMISSIONDIRECTION")
		{
			dir = get_hpv3d(ss);
		}
	}

	ion_source_ = new ParticleSource(type, voltage, max_emissions, center, dir);
}





void plotVectorField(char const * file_name)
{
	FILE *pipe = popen("gnuplot -persist 2>/dev/null", "w");
	fprintf(pipe, " set ylabel 'Y-AXIS'\n");
	fprintf(pipe, " set xlabel 'X-AXIS'\n");
	fprintf(pipe, " set zlabel 'Z-AXIS'\n");
	fprintf(pipe, " set yrange [ -.6 : .6 ]\n");
	fprintf(pipe, " set xrange [ -.6 : .6 ]\n");
	fprintf(pipe, " set zrange [ -.6 : .6 ]\n");

	std::stringstream ss("splot './");
	//"splot './magfile.tsv' with vectors filled head lw 2\n"
	ss << "splot '" << file_name << "' with vectors\n";
	fprintf(pipe, "%s", ss.str().c_str());
	pclose(pipe);
}

}
